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ABSTRACT 

We present results of a fully non-local model of convection for white dwarf envelopes. We 
show that this model is able to reproduce the results of numerical simulations for convective 
efficiencies ranging from very inefficient to moderately efficient; this agreement is made more 
impressive given that no closure parameters have been adjusted in going from the previously 
reported case of A-stars to the present case of white dwarfs; for comparison, in order to match 
the peak convective flux found in numerical simulations for both the white dwarf envelopes 
discussed in this paper and the A-star envelopes discussed in our previous work requires 
changing the mixing length parameter of commonly used local models by a factor of 4. We 
also examine in detail the overshooting at the base of the convection zone, both in terms of 
the convective flux and in terms of the velocity field: we find that the flux overshoots by 
^ 1.25///) and the velocity by ^ 2.5Hp. Due to the large amount of overshooting found at 
the base of the convection zone the new model predicts the mixed region of white dwarf 
envelopes to contain at least 10 times more mass than local mixing length theory (MLT) 
models having similar photospheric temperature structures. This result is consistent with the 
upper limit given by numerical simulations which predict an even larger amount of mass to 
be mixed by convective overshooting. Finally, we attempt to parametrise some of our results 
in terms of local MLT-based models, insofar as is possible given the limitations of MLT. 
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1 INTRODUCTION 

Understanding white dwarf stars is crucial for many areas of as- 
trophysics. First, their masses can be used to place constraints on 
mass loss in the post-Main Sequence phase and hence on the initial- 
final mass relation (Weidemann 2000). Second, their temperatures 
can be used to derive ages , either ind ividually, for clusters, or for 
the local Galactic disk (e.g. lwinget"et al. 1987; Wood 1992). Com- 
plementary to this, white dwarfs are observed to pulsate in spe- 
cific temperature ranges, and these pulsations all ow us to probe 
and constrain the interior structure of these stars jWingej |l998). 
Through their pulsations, we can use asteroseismology t o examine 
various physical processes such as nuclear reac tion ra tes iMetcalf j 
|2p03), chemical diffusion ( Montgomerv et al. 2001), crystallisa- 
tion tWinget et al^.., 1997; Montgornerv & Wingeliil99ft) . and neu- 
trino emission JO'Brienetalll998l) . 

In addition, these pulsations are most likely driven through 
their interaction with the surface convectio n zone in these stars 
(Brickhill 1991a; Wu 1997; Goldreich & w3ll999h . meaning that 
the onset of pulsations (in terms of T^a) is linked with the con- 
vection zone reaching a certain depth. In many cases the observed 



amplitudes imply that the depth and structure of these convection 
zones should vary appreciably (by a factor of several in mass) 
during a pulsational cycle, and this is thought to be the origin of 
the dominant nonlinearities in many of the observed lightcurves 
(Brickhill 1992b; Wu 2001). Given the simplification afforded by 
the separation of the convective turnover timescale 1 s) and the 
pulsation timescale (~ 100 s), it is possible to use the pulsations 
themsel ves to sample and co nstrain the convection zones of these 
stars ( Ising & Koester 2001, Montgomery, in preparation). Con- 
sequently, other than the well-studied convection zone of the Sun, 
white dwarfs may offer the best chance for testing theories of stellar 
convection. 

In this paper, we apply the Reynolds stress model formalism 
for convection (see Kupka & Montgomerv "2002") to envelope mod- 
els of both DA (hydrogen) and DB (helium) white dwarfs. As was 
the case for convection in A-stars, this problem is made easier nu- 
merically by the fact that the white dwarf models we consider have 
relatively thin surface convection zones, and hence shorter thermal 
relaxation timescales. As a consequence, we have treated models in 
which convection is not the dominant form of energy transport, i.e.. 
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Figure 1. (a) The fraction of the flux carried by convection for five DA white dwarf models with the indicated effective temperatures, where we have taken 
logT as our radial variable; logg = 8.0 and X = 1.00 (Z = 0.00) for all models. The cross on each curve near logT ~ 4.1 shows the location where r = 2/3 
for each model. As expected, convection becomes more dominant with decreasing T^ff. (b) The same as (a) but for the rms of the vertical component of the 
convective velocity. 



Fc/Ft < 0.5; in this regime, convection is much more sensitive to 
details in the modelling. 

In the sections which follow, we give a brief outline of the 
physics and the numerical procedure used to compute our envelope 
models. Results are then presented for models of DA and DB white 
dwarfs, as a function of Teff. Finally, where possible we compare 
our results to those of numerical simulations and to models used for 
fitting white dwarf spectra. We show that our calculations qualita- 
tively follow the simulations, especially in terms of the convective 
flux, and that the photospheric temperature structure of our models 
is similar to that of previous investigations. 



2 DESCRIPTION OF MODEL 

The convection mo del used here is an extension of the 
ICanuto & Dubovikovl (1998, hereafter CD98) model which re- 
quires the solution of five differential equations of first order in 
time and second order in space for the hydrodynamic moments K, 
9^, J = w6 = Fc/ipCp), W-, and e, and of an additional equation for 
the ti me evolution of T (cf. equations (l)-(5) and (8) in .Kupk3 
Il999h : here and throughout this paper, K is the turbulent kinetic 
energy, 6 and w are the temperature and velocity fluctuations, re- 
spectively, and e is the turbulent kinetic energy dissipation rate. 
This system is completed by an equation for the total pressure ('hy- 
drostatic equilibrium' including turbulent pressure, equation (7) in 
Kupka 1999b) and for the mass ('conservation of mass'). We solve 
this set of differential equations on an unequally spaced mass grid, 
with the zoning chosen so as to resolve the gradients in the var- 
ious quanti ties. The model equ ations were derived in a series o f 
papers bv iCanutd <199lll993h . CD98, and lCanuto et^ <200ll) 
Comp ressibility effects are taken into account following ICanutd 
ll993|). The adoption and solution of these equations for stellar 
envelopes is described in Kuoka & Monteomerv (2002), and it is 
this model which we apply here as well. The closure parameters 
for correlations such as 6dp' /dxj, i.e. between fluctuations of the 
temperature and the pressure gradient, which appear in the non- 
local model, have been taken over from the previous paper (see 
iKupk a & Montgomerv..200Z) . We note here that no mixing length 



is used in this model due to a separate evolution equation for the 
turbulent kinetic energy dissipation rate e. 

For our calculations, we have used a Prandtl number of 10"* 
as a typical value for the outer part of white dwarf envelopes; val- 
ues 2 orders of magnitude smaller than this produce essentially 
identical results. We note that values 2 orders of magnitude larger 
(i.e., 10"*) can alter our results for the hottest models at up to the 
~ 10% level. For the constitutive physics, we use the equation of 
state and opa city data from the OPAL project (Rogers et al. 1996; 
llglesias & Rogersi )1996). Since we are treating white dwarfs with 
pure or nearly pure surface layers, we have taken the metallicity 
to be zero (Z = 0.0), with X = \.Q for the DA (hydrogen spec- 
trum) white dwarfs and 7 = 1.0 for the DB (helium spectrum) white 
dwarfs. 



While the convection z ones we have treated here and in 
iKupka & Montgomery! (2002) are relatively thin surface convec- 
tion zones, we have solved the full equations for spherical geome- 
try. For upper (outer) boundary conditions, we fix the temperature, 
gravity, and stellar radius to be equal to those obtained from an en- 
velope model assuming local (MLT) convection, which itself has 
a given Tcff, logg, and R^,; thus, both the location of the model in 
the H-R diagram and its mass are specified (these may be taken 
either from a self-consistent stellar model or freely specified). We 
mention that the outer photospheres of these envelope models are 
completely radiative, so the use of MLT in these models leaves no 
direct imprint on our subsequent solutions, with the exception of 
the value derived for the stellar radius, which is slightly affected. 
At the lower boundary, we assume a constant input luminosity L* 
equal to the luminosity at the stellar surface. The complete system 
is then integrated in time (currently by a semi-implicit method) un- 
til a stationary, thermally relaxed state is found. The mass shells can 
be re-zoned to a different relative size to resolve, e.g. steep temper- 
ature gradients that may appear and/or disappear during conver- 
gence. The radiative envelope below the convection zone may then 
be obtained from a simple downward integration. 
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Table 1. Convection zone parameters of DA white dwarf models obtained witli the non-local mode l. The overs hooting (OV) is defined as the distance in 
pressure scale heights from the minimum of Fc/Fj to the point where |Fc/^t| ~ 10"^ (similar to Lt in Zahn 1991*). OV[mix] is the velocity overshooting; it is 
defined as the distance from the base of the formally unstable convective region (V ~ V^d) to the point where a linear extrapolation of the velocity vanishes 
(see Fig.|2). Finally, a^ff is the value of a for the given mixing length theory (either MLl or ML2) required to reproduce the peak flux in the convection zone. 



Teff logg log(Mcz/M*) (Fc/FT)max OV (vc)max (vc)t-=2/3 OV[mix] (Plurb/Ptol)max (/'turb/Plot)r=2/3 "eff 

(K) (in Hp) (kms-') (km.s-') (in Hp) (MLl) (ML2) 



12200 


8.00 


-14.72 


0.4402 


1.41 


6.43 


6.22 


2.75 


0.290 


0.242 


1.66 


0.65 


12400 


8.00 


-14.93 


0.2124 


1.37 


5.61 


4.84 


2.45 


0.172 


0.155 


1.57 


0.64 


12700 


8.00 


-15.07 


0.1020 


1.24 


5.30 


4.05 


2.22 


0.125 


0.106 


1.73 


0.72 


13000 


8.00 


-15.14 


0.0642 


1.21 


5.40 


4.22 


2.19 


0.125 


0.107 


2.04 


0.85 


13400 


8.00 


-15.22 


0.0346 


1.14 


5.50 


4.41 


2.15 


0.125 


0.106 


2.53 


1.06 


12700 


8.30 


-15.28 


0.3414 


1.29 


5.58 


4.98 


2.50 


0.184 


0.161 


1.54 


0.62 


12700 


7.70 


-14.70 


0.0473 


1.23 


5.88 


4.91 


2.29 


0.152 


0.141 


2.30 


0.96 



3 RESULTS 

In the following sections, we report results for DA (hydrogen) 
white dwarfs and DB (helium) white dwarfs. For most of our cal- 
culations we have used the canonical value of the white dwarf sur- 
face gravity of log^ = 8.0, which corresponds to a stellar mass of 
~ 0.6Mq, although we also present selected results for logg = 7.7 
and 8.3. 

The temperatures we have explored in these models range 
from high temperatures in which convection is very inefficient to 
cooler temperatures for which the convection becomes deeper al- 
though not yet adiabatic. Since convection in the DA's is driven by 
Hi partial ionisation and in the DB's mainly by He II partial ion- 
isation, the temperatures of the DB models are much higher; for 
our DA models, we have chosen temperatures in the range 12,200- 
13,400 K, and for the DB models 28,000-35,000 K. Observation- 
ally, these temperatures approximately correspond to the onset of 
pulsations in these stars, the so-called 'blue edge' o f the pulsational 
instab ility strip, which is ~ 12,500 K for the DA's {B ergeron et alJ 
Il995h and ~ 28,000 K for the DB's iBeauchamp et aljl99gl) . These 
temperatures can be explained either in terms of linear instabil- 
ity due to convective driving (Brickhill 1991a; Goldreich & Wu 
^999 ) or the traditional k-j driving mechanism ( Winget et alil982t 
rWinget.1982') . Thus, the nature of convection in these stars has a 
large impact on the properties of the pulsations. 



3.1 DA models 

For the DA's, we have considered temperatures between 12,200 K 
and 13,400 K; Fig.Qshows our central results. First, the models are 
all strongly convective in the photosphere, with the convective flux 
being a substantial fraction of the maximum value in the H I con- 
vection zone. Second, we see that the photospheric velocities are 
even larger, attaining values at least as large as 75% of their maxi- 
mum values within the convection zone. Far out in the photospheres 
of these models (the crosses indicate the point at which r = 2/3), 
we see that the models do become radiative, justifying our use of 
fully radiative outer boundary conditions. We note that the wiggles 
near \ogT ~ 4.35 and 4.07 in the velocity field are due to terms in 
the equations for third order moments, which represent non-local 
transport, and whose functional form depends on the sign of A'^ 
(the square of the Brunt- Vaisala frequency); a smoother model for 



these terms would remove this feature from the velocities.' The 
actual velocity distribution in transition regions is expected to be 
smooth, and is found to be so in numerical simulations. We also 
have to emphasise that the log T scale does not properly resolve the 
upper and mid photosphere, hence the apparently very steep drop 
of velocity at the surface. Comparisons with the convective flux 
displayed in Fig.^ and the kinetic energy flux shown in Fig.|3|be- 
low indicate what a plot as a function of logarithmic optical depth 
(logTross) reveals more clearly: a steady and by no means abrupt 
decay towards zero. Moving the upper boundary further outwards 
by two orders of magnitude on the optical depth scale would not 
be noticed in Fig.0), as the photospheric temperature has already 
reached a constant value there due to the underlying approximation 
of grey radiative transfer. 

In Table □ we give a summary of our results for these DA 
models. For each model, we give the maximum value of the con- 
vective flux, (Fc / FT)m-dx, the overshooting in pressure scale heights 
of the flux, OV, and the velocity, OV[mix], and the maximum and 
photospheric values of the vertical component of the convective ve- 
locity and the turbulent pressure, vc and ptuit/ptot, respectively. We 
also give log(Mcz/A^*), where Mcz is the total mass of the convec- 
tion zone, defined as the region which is mixed (including velocity 
overshooting). Finally, we list values of Oeff, which is the value 
of Q which a given MLT model requires in order to reproduce the 
maximum convective flux. The columns labelled 'MLl' and 'ML2' 
denote the values obt aine d using two different versio ns of MLT, by 
lB6hm-Vitens3 <1958^ and'Bohm & Cassinellil jl97lh . respectively. 
These versions of MLT will be discussed further in Section|4| 

Over the temperature range of the models, we see that for MLl 
convection that a^f lies roughly in the range 1.5-2.5; clearly, a lo- 
cal MLT using a single value of a would be unable to reproduce 
these results. At any rate, this range of values is consistent with 
those found by Ludwig et al. ( 1994) in comparison with numerical 
simulations, and is also consistent with the model atmosphere fits 
of Koester et alj tl994') . In Section|4| we make a more systematic 
comparison with MLT models (in terms of both MLl and ML2), 
taking into account the photospheric temperature structure of the 
models. 



For instance, the non- Gaussian closures recently suggested by 
iGrvanik & Hartmand I200I should not produce such a feature. Their new 
model promises a more realistic approach to third and fourth order moments 
in Reynolds stress models, but requires the solution of additional differen- 
tial equations. 
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Figure 2. A comparison of the convective fluxes (upper panel) and veloci- 
ties (lower panel) of the non-local model and of an MLT model (/ = 1.73//p), 
for Teff = 12,700 K. We see that the MLT model cannot simultaneously 
match the height and width of the convective flux. In addition, the veloc- 
ity field of the MLT model does not extend as far out into either the pho- 
tosphere or downward into the envelope as that of the non-local model. 
Along the top axis we have indicated the depth in terms of the envelope 
mass log(l -Mr/Mt), so that we can see the relative amounts of mass in the 
formally unstable and overshooting regions. 



In Fig.|2| we make a detailed comparison between the flux and 
velocity structure of our Teff = 12,700 K envelope model with that 
obtained using the standard MLT prescription for convection. There 
are at least four aspects of our convection zone solutions which we 
might wish to compare to MLT; the peak flux and its overshooting 
depth (i.e. the convectively mixed region), and the peak velocity 
and its overshooting depth. Given the inadequacies of a local MLT, 
however, it is only possible to fit one of these quantities at a time. 
Thus, we have adjusted the mixing-length parameter a so that the 
two models have the same maximum flux; this is achieved for a 
value of Q = 1 STiHp. In terms of the convective flux, we see that the 
non-local model predicts a wider convective region with overshoot- 
ing extending down to logT ~ 4.53, as compared to the local MLT 
model, which has its base at logT ~ 4.3. Significantly, the non- 
local model also predicts a much larger value for the photospheric 
flux. In terms of the convective velocities, the differences are even 
larger. The non-local model predicts velocities 50% larger than the 
local model, and due to overshooting, these velocities also extend 
much deeper, down to log T ~ 4.6, as compared to log T ~ 4.33 for 
the local model. Thus, we find that the mass of material 'stirred' 
by the convection zone is an order of magnitude larger than that 
which is formally unstable according to a local stability criterion. 
Among other things, this can have consequences for the diffu sion of 
eleme nts in white dwarf envelop es, as was pointed out bv iFrevtad 
il995h and lFrevtag et alj i 1996ft . who, in addition, found an even 
larger amount of material stirred below the convection zone in their 
numerical simulations of DA envelopes (see also Section|5j. 





- 1 1 1 1 1 








r 

















X " 








^ - 


0.005 


\ 

\ 








- 

- 


-0.01 


- 

- 








- 

- 




: 12200 


A 






- 


0.015 


- 12400 


K \ 










; 12700 


K \ 










: 13000 


K 








-0.02 


- 13400 


K 










■ 1 , , , 1 












4.8 4.6 




4.4 4.2 




4 



log T 

Figure 3. The kinetic energy flux as a function of log T, for the five different 
models shown in Fig.0 As expected, the cooler models have larger fluxes. 
Most significantly, however, these numbers show that |Fi[i„| is essentially 
negligible for the models we have examined. 



3.2 Kinetic energy flux and non-locality 

In Fig.|3|we plot the kinetic energy flux for the five different mod- 
els shown in Fig. Q Besides the fact that the cooler models have 
larger fluxes, which is to be expected, we see from the magnitudes 
of these fluxes that Fyn is essentially negligible for the models we 
have examined. Taking this result at face value, we are thus in a 
different regime from that of the Sun, in which |Fki n/FT| may be 
as lar ge as 20 per cent (cf. IStein & Nordlund.l993:lKim & ChanI 
Il998t) . Although a positive Fun in the photosphere indicates that 
the skewness of spectral lines (and their NLTE cores) produced 
in this region sho uld also be positive (CD98, cf. the discussion in 
iKupka & Montgom erv 2002), the small magnitude of Fkin means 
that this effect will be small and probably difficult to measure. In 
the future, it will hopefully be possible to make observational tests 
of these predictions. 

Although the kinetic energy flux in these models appears to 
be small, this transport of kinetic energy may lead to an equilib- 
rium state having large velocity fields outside of the formally con- 
vective regions, provided that the local dissipation rate of kinetic 
energy is small enough. For instance, in Fig. |2|we see that both 
below (logT > 4.4) and above (logT < 4.05) the formally con- 
vective region there is significant overshooting of the velocities. 
This can be understood from Fig.|3]since Fki„ is positive at the top 
of the formally convective region (logT ^ 4.1), indicating that ki- 
netic energy is being transported outward, while Fkin is negative at 
the bottom of the formally convective region (log T ~ 4.4), indicat- 
ing that kinetic energy is being transported inward. This velocity 
overshooting is an essential feature of these models, and is an ef- 
fect which cannot be captured using local MLT-type models (see 
Fig.|2|above). 

Moreover, even in the formally unstable part of the convection 
zone, where Fki„ is small in magnitude compared to Fco„v (around 
logT < 4.2, i.e. where the largest velocities occur), non-local ef- 
fects remain important. The reason for this result is that in the gov- 
erning equation for the velocity field in the non-local model, the 
third order moment directly related to Fki,,, q-w = 2Fkin//5, appears 
under a divergence (eq. 1 in|KuDka 1999, eqs. 19a, 36a, and 51a in 
CD98). Thus, the rapid variation of q-w can cause its divergence 
to be large even if the third order moments themselves are small 
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(a) log T (b) log T 

Figure 4. (a) The fraction of the flux carried by convection for six DB white dwarf models with the indicated effective temperatures, as a function of logT; 
logg = 8.0 and Y = 1.00 (Z = 0.00) for all models. The crosses in the range of log T ~ 4.45^.55 show the respective locations of the photospheres (r = 2/3) 
of the different models. As expected, convection becomes more efficient with decreasing T^ff. (b) The same as (a) but for the rms of the vertical component of 
the convective velocities. 



Table 2. The same as TablelTI but for DB white dwarf models. 


^cff 




\og(Mcz|M^,) 


(Fc/FT)max 


OV 


(Vc)max 


(vc)t=2/3 


OV[mix] 


iPtuxbl /'tot)max 


(/'lLrb/Aot)r=2/3 




(K) 








(inffp) 


(km S-' ) 


(kms-') 


(in//,,) 






(MLl) 


(ML2) 


28000 


8.00 


-13.28 


0.4030 


1.33 


5.67 


4.67 


2.53 


0.168 


0.166 


1.14 


0.45 


29000 


8.00 


-13.42 


0.2577 


1.32 


5.46 


2.62 


2.36 


0.102 


0.057 


1.07 


0.43 


30000 


8.00 


-13.48 


0.2086 


1.28 


5.72 


2.24 


2.36 


0.111 


0.040 


1.12 


0.45 


31000 


8.00 


-13.51 


0.1683 


1.29 


5.99 


2.04 


2.36 


0.120 


0.032 


1.16 


0.48 


33000 


8.00 


-13.55 


0.1136 


1.31 


6.73 


1.93 


2.47 


0.147 


0.027 


1.28 


0.53 


35000 


8.00 


-13.58 


0.0751 


1.34 


7.60 


1.71 


2.56 


0.177 


0.020 


1.42 


0.59 


33000 


8.30 


-13.96 


0.1696 


1.31 


6.03 


2.18 


2.41 


0.120 


0.034 


1.20 


0.49 


33000 


7.70 


-13.12 


0.0736 


1.30 


7.69 


1.51 


2.56 


0.182 


0.016 


1.39 


0.58 



(e.g. as compared to the convective flux). This divergence couples 
the velocity field with the vertical structure of the model in an in- 
teresting manner: fluid elements having differing degrees of par- 
tial ionisation are transported to and from neighbouring layers. The 
ionisation and recombination work and the ionisation energy con- 
tained in the gas change the nature of the flow. They hence play an 
important role in the non-local nature of convection driven by ioni- 
sation zones (for the related case of numerical simulations of solar 
granulation see Rast et al. 1993). 

As a consequence, for models of both 'early' white dwarfs 
and early A-stars we find that their convection zones can be in- 
efficient in the sense of having small convective and kinetic en- 
ergy fluxes while still having quite large velocity fields, as seen in 
Fig.Q and that such large velocities can have significant observa- 
tional and theoretical consequences. Fo r A-stars, the studies of line 
profiles discussed in lLandstreell i 19981) provide a strong observa- 
tional indication for large velocity fields with a large asymmetry 
between up- and downflows in these objects. Although such ob- 
servational indicators are not yet available for white dwarfs, the 
similarity of non-local models for both A-stars and white dwarfs, 
which is suppo rted independently by the numerical simulations of 
iFrevtad i 19951 see also Section|5), suggests that white dwarf en- 
velopes can have significant convective velocities despite having a 
nearly radiative temperature structure. Thus, the non-local effects. 



far from being unimportant, manifest themselves in relatively large 
convective velocities, which in turn can affect the formation and 
shape of spectral lines in the atmosphere, the diffusion/gravitational 
settling of elements in the envelope, as well as the way in which any 
pulsations present in the star interact with the convection zone. 



3.3 DB models 

For the DB's, we have calculated models with temperatures be- 
tween 28,000 K and 35,000 K. The analogous results to the DA case 
are given in Fig.|4|and in Table |2| Besides the general trend that 
convection becomes more efficient as the temperature decreases, 
we see that these results are actually quite similar to what we 
found for the DA's. For instance, the maximum convective veloc- 
ities of both are in the range 5-6 km/s, their flux overshooting is 
~ 1.25 Hp, and their velocity overshooting is ^ 2.5 Hp. One dif- 
ference between the DA and DB models is the remarkable increase 
in the photospheric velocity seen in Fig05 for the coolest model 
(28,000 K). This is caused by additional driving provided by the 
partial ionisation of He I. At such low optical depths (t ~ 2/3) this 
ionisation zone cannot alter the energy transport, but interestingly, 
at least within the context of the present, non-local model, it is able 
to leave a very clear feature in the photospheric surface velocity 
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Figure 5. Upper panel: The temperature structure of the photosphere of the 
r(.ff=12,200 K, logg=8.0 model as a function of optical depth of the non- 
local model (solid curve) compared to that of local MLT models (MLl) 
having a = 1.5 (dotted curve) and a = 3.0 (dashed curve). Lower panel: 
The residuals from a least squares fit as a function of a over this range. The 
best fit is achieved for a =2.23; for ML2, we find a =0.89. 

field. For higher effective temperatures, this ionisation zone moves 
to even smaller optical depths and its influence on convection van- 
ishes completely. 

The other principal differences between the DA and DB mod- 
els turn out to be in the depth of their convection zones and in their 
convective efficiencies, as parametrised by Qcffi the DB's have sys- 
tematically lower values of Ocff, and for a given convective flux the 
DB's have convection zones which are an order of magnitude more 
massive than those of the DA's. As it turns out, these two results 
are consistent with each other. Since the He II convection zone is 
deeper than the H I zone, the MLT eddies contain denser material 
which has a higher heat capacity and is optically thick, allowing 
them to transport the same fraction of the total flux with a smaller 
value of Q, i.e., more efficiently. 



4 COMPARISON WITH MLT 

There are many different ways in which a comparison between the 
Reynolds stress model and local MLT can be made. In Section^ 
we chose to match the maximum of the convective flux in the two 
models. Since the maxima occur at large optical depths, we are in 
some sense matching the models in a 'deep' region. Alternatively, 
we could choose to do the comparison in the photosphere, since this 
has relevance for model atmosphere fits. In this section, we com- 



pare the temperature structure of the photospheres of the local and 
non-local models; a comparison of the velocity fields is not infor- 
mative since the local models produce convection zones which are 
always too narrow for a given convective flux. In addition, the ve- 
locity fields used in model atmosphere fits are not usually derived 
from an underlying convection model but are treated as a free pa- 
rameter ('microturbulence'), whereas the temperature profiles used 
in the fits are computed using a model of the convection. 

Since the temperature structure is most dependent on the con- 
vective prescription for models with the largest photospheric con- 
vective fluxes, we choose our lowest temperature models for the 
comparisons. In Fig.|5]we show the results of such a fitting proce- 
dure for our rcff=12,200 K, logg=8.0 DA model. In the upper panel 
we plot the temperature as a function of optical depth. We see that 
the non-local model (solid curve) is bracketed by the local MLT 
models (MLl) having a = 1.5 (dotted curve) and a = 3.0 (dashed 
curve). In the lower panel, we plot the results of a least squares fit 
of the MLT profiles to the non-local profile, as a function of a. The 
best fit is achieved for a value of a = 2.23. This is in broad asree- 
ment with Koester et al. ( 1994^, who find that a = 2.0 yields model 
spectra w hich are in reasonable agreement with the observations. In 
addition, (Ludwig_et al. ( 1994) found that the temperature structure 
of their reff= 12,600 K hydrodynamic model could be approximated 
with values of a between 1 and 2. For our reff=12,600 K model, 
we find it to be bracketed by a = 1 .0 and 2.5 models, with a best fit 
value of a = 1.86, which is consistent with their result. 

Many of the current model atmosphere fits for white dwarfs 
are based on a version of MLT i n which radiative energy losses 
are somewhat suppressed (Bohm & Cassinellilll97lh ; this version, 
which is more efficient than MLl for a given val ue of a, is ofte n 
referred to as ML2 (e.g. iTassoul e t al. 1990; Lud wig et alJll999l) . 
For our 7'eff=12,200 K model, we find that it is best fit with a value 
of a = 0.89. For our reff=12,600 K model, we find that the best fit 
value of a is 0.73. This is in good agreement with Bergeron et al] 
jl995). who found that their optical spectra of ZZ Ceti's could be 
well fit by ML2/a = 1 .0 models, but that best fits to both the UV 
and optical spectra were obtained with ML2/a = 0.6, which is the 
value they adopted for their subsequent analysis. 

For our DB models, we are unable to do such fits of the pho- 
tospheric convective efficiency since the models have so little con- 
vection in the photosphere. Thus, the only results comparing our 
models with local MLT models are the ones already given in Ta- 
ble|2| 



5 COMPARISON WITH NUMERICAL SIMULATIONS 

In Table|21 we compa re the res ults of our DA models to those from 
the 2D simulations o f lFrevtaja995.) . This comparison is not com- 
pletely fair due to the fact that different equations of state have been 
used. Even so, we see that the agreement in terms of the maximum 
convective fluxes is actually quite good, although the photospheric 
fluxes do not agree as well. In addition, both the maximum and 
photospheric values of our vertical convective velocities (vc) are 
systematically higher than his by a sizeable margin. At the present, 
we do not know if this signals an incompleteness in the Reynolds 
stress m odel or an inade quacy in the 2D calculations, or both. For 
instance. I Asplund et alJ )2000) found that while 2D simulations in 
the Sun did reasonably well at reproducing the temperature struc- 
ture of the 3D simulations (and therefore the convective fluxes), the 
2D simulations systematically underestimated the magnitude of the 
convective velocities by 10-20%. In addition, the viscosity used for 
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Table 3. Comparison of Reynolds stress model results to the numerical simulations o f lFrevtad IT99I : all models 
have log g=8.0. We have used vq and \\ to denote the vertical and horizontal rms convective velocities, respectively. 
Columns labelled RS are the Reynolds stress results, and columns labelled F are the results of Freytag's numerical 
simulations. 
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the 2D the simulations may also lead to a reduction of the veloci- 
ties. Conversely, it is possible that some of the closures in our im- 
plementation of the Reynolds stres s model may be responsible for 
this discrepancy (e.g. see Fig. 5 of lKupkjEool . Finally, perhaps 
surprisingly, we note that the discrepancy between our models is 
smaller for the horizontal component of the convective velocity, Vh, 
than it is for the vertical component. 

An interesting question con cerns how the vel ocity field de- 
cays beneath the convection zone . iFrevtad \ 1 99^) and lFrevtag et alj 
^996) claim both numerical and theoretical evidence for an expo- 
nential decay with depth of the turbulent velocity field; their posi- 
tion is given some support by Ludwig ( 2003) for the case of over- 
shooting into stellar atmospheres, although he does not find this to 
be the case for all the models examined. Our equations, on the other 
hand, lead both analytically and numerically to an approximately 
linear decay of the velocity field with depth. We claim that the cor- 
rect behaviour at the base of a convection zone is not yet known, 
since the way in which one filters out the travelling waves from the 
truly convective fluid motions has an important effect on the ve- 
locity field inferred from the hydrodynamical simulations (Ludwig 
I2OO3I) . Waves, on the other hand, are much less efficient in mix- 
ing a fluid than a network of drafts and plumes associated with a 
'genuine' convective velocity field. The resolution of this problem 
has an important bearing on our understanding of diffusion in white 
dwarf envelopes, since the depth to which convective fluid motions 
penetrate can greatly affect the diffusion of chemical elements. 

Finally, we mention a technical point concerning the Reynolds 
stress mod el, which we examine in det ail in an appendix. The for- 
malism of Canut o & Dubovikovl (l99^, which we employ, implic- 
itly assumes for the purpose of statistical averaging that the equa- 
tion of state is that of an ideal gas. This means that thermodynamic 
quantities are held constant during averaging, i.e., cpw9 is replaced 
by cpwO} This is not strictly valid since convection zones usually 
coincide with regions which are partially ionised, so that cp (and 
other thermodynamic quantities) are not constants but functions of 
temperature and density. In appendix]^ we examine the effect of 
this approximation and show that it leads to errors in the convec- 
tive flux no larger than 10-20% for the calculations we have done. 
Since this uncertainty is less than that introduced by other aspects 
of the modelling, such as d ifferent p rescriptions for the third order 
moments (see Fig. 5 o flKuDkal2003l for the case of A-stars), we are 
justified in neglecting this effect in the present analysis. In addition, 

^ We note that in our subsequent implementation of the Reynolds stress 
equations that we do use the actual values of cp and other thermodynamic 
quantities computed using the OR^L equation of state. 



we show that in the future the correction terms arising from relax- 
ing this assumption can be naturally included within the Reynolds 
stress formalism. 



6 COMPARISON WITH PULSATION DATA 

For pulsations in the DA and DB stars, we are in the fortunate 
regime in which the convective turnover time is quite short (~ 1 s) 
compared to the periods of the observed modes (~100 s). This 
means that the convection zone responds to the pulsations in a 
quasi-static manner, so the only input needed for the pulsation cal- 
culations is a sequence of static convection zones computed for dif- 
ferent effective temperatures. Thus, true time-dependent solutions 
for the convective region are not needed, which simplifies the mod- 
elling considerably. 

It has recently been shown terickhilll Il991albt 
iGoldreich &"wi] Il999ll that the convection zone plays a vital 
role in the driving of pulsations in the DAV stars, and this is likely 
to be true of the DBV stars as well. Thus, t he observed blue edge 
and r ed edge of the instability strip (see iBergeron et aljlT995l 
l2003t) should contain information concerning the dependence of 
the convective efficiency on Tcff and logg. 

In addition, the observed nonlinearities in the lightcurves of 
the DAV's are also believed to b e due to the interaction of the con - 
vection zone with the pulsations l'Wu'200il llsing & Koester*2 00lh . 
so it should again be possible to use the observations to constrain 
different models of convection; this work is presently under way 
(Montgomery, in preparation). Thus, the pulsating white dwarfs 
(both DAV's and DBV's) offer a great deal of promise for learn- 
ing about the physics of convection. 



7 CONCLUSIONS 

Using a fully non-local model of convection together with a real- 
istic equation of state and opacities, we have calculated envelope 
models for stellar parameters appropriate for DA and DB white 
dwarfs. We find good agreement between our models and those 
obtained through fitting white dwarf spectra, as well as good agree- 
ment with the results of hydrodynamic simulations. 

First, the maximum convective fluxes in our DA models com- 
pare reas onably well w ith those found in 2D hydrodynamical sim- 
ulations iFrevta di l995h . This result appears more impressive when 
taking into acc ount that a similar agreem ent was found for the case 
of A-stars I Kupka & Montgomervll20()31 while MLT type models 
such as MLl require a change of the scale length parameter a 
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by a factor of 4 (see already iFrevtadfTQQSh . On the other hand, 
the convective velocities (both vertical and horizontal) differ from 
each other by a fairly significant amount. At present we do not 
know whether this signals an incompleteness in the Reynolds stress 
model or a limitation of the 2D simulations, or both. However, 
while our rms velocities and those of Freytag (1995) are found to 
differ by up to a factor of two in the convection zone, we empha- 
sise the completely different results obtained when using a local 
convection model such as MLl or ML2. Among other shortcom- 
ings, the latter underestimate the size of the convectively mixed 
region by a factor of at least 10 in terms of mass re lative to our new 
results and even more when compared to those of Frevtag ( 1995). 
Thus, while the predicted rms velocities at least qualitatively agree 
when comparing the Reynolds stress formalism applied in this pa- 
per with 2D numerical simulations, results from local convection 
models such as MLT are fundamentally different from the two. 

Second, given the widespread use of MLT in stellar and at- 
mospheric modelling, we have compared aspects of our models to 
MLT models. We find for our DA models that their photospheric 
temperature structure is best approximated by local MLT models 
(ML2) having a between ~ 0. 7 and ~ 0.9. This is in agreement with 
the results of iBergeron et alj Jl995h . who found that their model 
spectra matched both the optical and UV data best for ML2/a=0.6; 
for fitting the optical data alone, ML2/a=1.0 also provided a good 
fit to their observations. Going beyond such fits, our non-local 
model offers the promise of a unified model for self-consistently 
computing the velocity field ('microturbulence') which is also re- 
quired for such model atmosphere fits, something which local mod- 
els fail to do. 

Since the onset of convection goes hand in hand with the on- 
set of pulsation in both the DA's and DB's, we have the opportu- 
nity to use asteroseismology to sample their interior structure. First 
of all, the dependence of the blue edge of the instability strip on 
Teff and logg very likely depends on how the convective efficiency 
varies with these parameters, so the observed blue edge can be used 
to constrain theories of convection. In addition, recent calculations 
have shown that the convection zones in these stars may b e the main 
source of the observed non-linearities in their lightcurves te rickhilll 
[19923; Wu 2001; Ising & Koester 2001), so it may be possible to 
use the pulsations themselves to probe the structure of the convec- 
tion zone. In particular, from observations of a given pulsating star 
it may be possible to infer the depth of the convection zone as a 
function of the instantaneous effective temperature (Montgomery, 
in preparation). This would provide important data for the current 
convection models. 

From an astrophysical standpoint, white dwarf stars are cru- 
cial for our understanding of the initial-final mass relation for 
stars as well as providing an independent method for determin- 
ing stellar and Galactic ages. In addition, through aster oseismol- 
ogy, t hey can serve as test b eds for nuclear reaction ra tes iMetcalf j 
|2003"), chemical diffusion ^Montgomerv et al."200l"), crystallisa- 
tion ( Winget et al. 1997; Montgomery & Winget 1999), and neu- 
trino emission I O'Brien et a l. 1998). In order to make maximum 
use of these inferences, however, model atmospheres need to be 
applied to the observations of individual white dwarfs in order to 
derive the stellar parameters (M*, Tcff); depending upon the pre- 
scription which is adopted for conv ection, these parameters can 
vary significantly (e.g. IBergeron et alll995h . Thus, understanding 
convection in the context of these stars is vital. 

Except for a possible application to RR Lyr stars and 
Cepheids, the results reported here and in lKupka & Monteomeryj 
j2ooa) represent the limits of what we can treat using our present 



semi-implicit solver. We are currently in the process of developing 
a fully-implicit solver, which will allow us to treat thicker convec- 
tion zones having larger thermal timescales, and, hence, lower ef- 
fective temperatures. In addition to allowing us to compute cooler 
models with deeper convection zones for the DA and DB white 
dwarfs, and the A and F-stars, we would eventually like to study 
the convection zone of the Sun. This would be an excellent test 
of the model since a wealth of high-quality data already exists for 
the solar convection zone. On the same basis, calculations of mix- 
ing and overshooting in stellar convective cores may be an even 
more important application. While the quality of data for probing 
the effects of convective core overshooting cannot match some of 
the solar observational results, there are no extra complications as 
introduced by the large temperature fluctuations and non-grey ra- 
diative transfer near the solar surface, which are unaccounted for 
by the present model but are expected to alter basic quantities such 
as the derived depth of the solar convection zone. The current accu- 
racy of evolution models of intermediate and high mass stars is still 
most severely limited by the coarse treatment of convective heat 
transfer and mixing. At least part of the solution to this problem 
might be achieved by using a model similar to the one applied in 
this paper. 
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APPENDIX A: ESTIMATE OF ERRORS ARISING FROM 
THE IDEAL GAS ASSUMPTION 

In this appendix, we give an assessment of the uncertainties intro- 
duced by the ideal gas assumption in the Reynolds stress model 
formalism of ICanuto & Dubovikov C1998,) we mentioned in Sec- 
tion|3 Essentially, this has meant that the specific heats, cp and cv, 
were treated as constants when computing ensemble averages. To 
partially correct for this, we have used a realistic equation of state 
to compute the values of cp and cv as functions of p and T where 
they occur in the relevant equations. However, we are necessarily 
missing any cross terms (moments) that may have arisen when they 
were averaged against other fluctuating quantities to derive the mo- 
ment equations themselves. 

To provide an estimate of the errors which this may introduce, 
we consider th e convective fl ux, which is actually the flux of en- 
thalpy, h. Rrom lCanutd i 1997I) . equation (19b), we have 

Ff = phuj, 

which, in the absence of a mean flow, we approximate as 



Ff = ph'u), 



(Al) 



in the Boussinesq sense. Here, Uj is the j'-th component of the tur- 
bulent velocity, and the overbar denotes an ensemble average. In 



the present context, we do not consider density-weighted averages, 
which are more convenient when dealing with a fully compress- 
ible flow. In what follows, all variables have been separated into 
an average and a fluctuating component with respect to ensemble 
averages. Ror example, if X represents a given fluctuating quantity 
(e.g. temperature, velocity, etc.), we write 



x = x+x', 

where 



(A2) 



X' = 0. (A3) 

In the work of Canuto fcanuto' 1997: ' Canuto & Dubovikov! 

1 1998ft the assumption has been made that the gas is ideal, specifi- 
cally that h = cpT, with cp being constant. Here, we consider h to be 
a function of T and P, i.e., h = h(T,P). If we assume that T' <C T, 
P' < P, we have 



h = h(T + T',P + P') 



= h{T,P) + hTT' +hpP' + -hTTT''' 



-0(T'p',P"), 



(A4) 

(A5) 



where 

hr 
hp 

Htt 
and 



Xp 



ah \ 
arjp 
Ml] 
dpJr 

d-h ' 



91np 
d\nT 



(1-vt)/p 

dT 



P''T 
T 



and each of the quantities hT,hp,hTT,Cp,h'T is evaluated at T.P.'p 
both here and in what follows. Since in the ideal gas case the first 
order temperature term is already present, and since we wish to go 
one order beyond this in the temperature and pressure fluctuations, 
we take the expansion to second order in temperature and first order 
in pressure. As we will see later, the first-order pressure corrections 
and the second-order temperature corrections are indeed of similar 
magnitude. 

Taking the ensemble average of equation <A5> . we find that 
h = h(T,P)+^hTTT^. 

Using h' = h- h, we can substitute the above result into equa- 
tion <AH for the convective flux. Writing w and 6 for the turbulent 
velocity and temperature perturbations, (i.e., 9 = T' , w = uj), we 
find that 

I 



Fc = phrwO + phpwP' + -phpT w{9^-9^) + ... 

Since hj = cp, the first term is the zeroth order term which we have 
been calling the convective flux, and the other terms are the correc- 
tions to it. Let us denote these terms by 



Fo 
Fi 



p Cp w9 
JihpwP' 



phTT w(9^-9-) 
J) Htt w9^ . 



(A6) 



The quantity wP' is not directly calculated in the present version 
of the code, al though it may be app roximated by closures such as 
those found in ICanutdiT99all997t) . Ror our present purposes, we 
use the relation 
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p' = 



dT , 



T = 



dP_ 

df 



to recast wP' in terms of w6, which would be exact in the incom- 
pressible limit. 

In Fig. lAll we examine the relative sizes of these terms for 
two different white dwarf models: a 12,200 K DA model (a) and a 
28,000 K DB model (b). In the upper panels we show the rms frac- 
tional temperature fluctuations as taken from our converged non- 
local models; we note that these fluctuations are everywhere less 
than 8%, justifying our treatment of them as small perturbations. 
In the lower panel, we plot the leading order flux, Fo (long-dashed 
curve), the first order pressure correction to it, Fi (dotted curve), the 
second order temperature correction, F2 (short dashed curve), and 
the sum of all of these. Ft (solid curve). We see that while these 
correction terms are non-negligible, they would lead to fractional 
corrections to Fo of only about 15%. In addition, these corrections 
appear to be large only in the central part of the convection zone; 
outside of this region, the convective (enthalpy) flux is essentially 
given by Fo. 

Finally, in Fig. IA2I we make a similar plot for an A- 
star model (reff=7200 K, log .g = 4.4) from our previous paper 
teupka 81 Montgomervl2o6^ . In contrast to the white dwarf case, 
we see that there are two different convectively unstable regions, 
corresponding to He II ionisation (log T ~ 4.7) and H I ionisation 
(logF ~ 4.1). In the deeper He II zone, we see that these correction 
terms are negligible. In the H I zone, these terms are larger, driven 
by the fairly large temperature fluctuations of up to ~ 18%. We see 
that the main effect is to shift the maximum of the convective flux 
slightly outward and to decrease its magnitude; even so, this de- 
crease in the maximum convective flux is only about 10%. As in 
the white dwarf case, we see that these corrections are only signifi- 
cant within the formal convection zone and not in the overshooting 
regions. 

We note that for both Figs. lAll and lA2l the models have not 



been reiterated, i.e., the converged model does not include the cor- 
rections to Fo. However, as the difference between Fo and Ft is 
no more than 15% and as there is a negative feedback, the figures 
are sufficient to demonstrate the size of the effect. Self-consistent 
models for the cases shown in Fig. lAll which are based on the full 
equation <A6> for Ft are expected to have a larger convective driv- 
ing and thus a maximum in Ft larger than the one plotted, though 
still smaller than Fo. The case shown in Fig. lA2l is not quite so sim- 
ple, but a self-consistent solution can be expected to yield a flux 
which lies between Fr and Fo, regardless of which one is larger. 

The additional correction terms which we have derived above 
(i.e., those in eauation lA6> pose no problem for the Reynolds stress 
approach since these new terms involve moments which are already 
calculated in the current formalism (wd,wd^ in the current paper, 
and eventually wP'). Thus, even if it turns out that these terms are 
more important for thicker/cooler convection zones, it should not 
prevent us from extending the Reynolds stress approach further into 
this regime. 

Finally, we note that the magnitude of the uncertainties in 
equation <A5t is of order ~ 10%, which is actually less than that 
introduced by other aspects of the modelling, such as differen t pre- 
scriptions for the third order moments (see Fig. 5 of ,Kupkai2003l 
for the case of A-stars). Thus, in the present analysis we are justi- 
fied in neglecting this effect. 
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White dwarf envelopes: further results of a non-local model of convection 
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Figure A2. Same as Fig. lAll but for an A-stai' model having Teff=7200 K 
and logg = 4.4. 
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